library(haven)
library(ggplot2)
library(tidyr)
#library(cshapes)
library(dplyr)
library(sqldf)
library(readxl)
library(scales)

setwd("C:/Users/jonves/Dropbox/Papers/SSP-Instability-Opinion/replication_data/")


gini <- function (x, weights = rep(1, length = length(x))) {
  ox <- order(x)
  x <- x[ox]
  weights <- weights[ox]/sum(weights)
  p <- cumsum(weights)
  nu <- cumsum(weights * x)
  n <- length(nu)
  nu <- nu/nu[n]
  gini <- sum(nu[-1] * p[-n]) - sum(nu[-n] * p[-1])
  return(gini)
}


ssp <- readRDS("data/ssp_long.rds")

models <- c("OECD Env-Growth", "IIASA GDP", "IIASA-WiC POP")
scenarios <- c("SSP1", "SSP2", "SSP3", "SSP4", "SSP5")
years <- c(2010, 2100)

variables <- c("region", "year", "scenario", "model", "value")

ssp <- ssp[which(ssp$model %in% models 
                 & ssp$scenario %in% scenarios
                 & ssp$year %in% years), variables]

df <- spread(ssp, model, value)


df$OECD <- df$`OECD Env-Growth`/df$`IIASA-WiC POP`*1000
df$IIASA <- df$`IIASA GDP`/df$`IIASA-WiC POP`*1000
df <- df[c("region", "year", "scenario", "OECD", "IIASA")]

df <- gather(df, model, gdpcap, -region, -year, -scenario)
df$model <- paste(df$model, df$scenario)
df <- df[c("region", "year", "model", "gdpcap")]
df <- spread(df, year, gdpcap)

names(df) <- c("region", "model", "gdpcap_2010", "gdpcap_2100")

ggplot(df, aes(x=log(gdpcap_2010), y=log(gdpcap_2100/gdpcap_2010))) + geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap(~model, nrow=2, ncol=5)

ssp32 <- readRDS("data/ssp_32_long.rds")
ssp32 <- ssp32[which(ssp32$variable %in% c("GDP|PPP", "Population") & ssp32$model == "PIK GDP-32"),]
ssp32$unit <- NULL
ssp32 <- spread(ssp32, variable, value)
ssp32$gdpcap <- ssp32$`GDP|PPP`/ssp32$Population*1000
ssp32$model <- paste("PIK", ssp32$scenario)
ssp32 <- ssp32[which(ssp32$year %in% c(2010, 2100)), c("region", "year", "model", "gdpcap")]
ssp32 <- spread(ssp32, year, gdpcap)
names(ssp32) <- c("region", "model", "gdpcap_2010", "gdpcap_2100")

df <- rbind(df, ssp32)

df$y <- log(df$gdpcap_2100/df$gdpcap_2010)
df$x <- log(df$gdpcap_2010)


# I know b coefficients are negative! Therefore, I transform to abs(coef), and use - instead of + in equation!
lm_eqn <- function(y, x){
  m <- lm(y ~ x)
  eq <- substitute(italic(y) == a - b * italic(x), 
                   list(a = as.numeric(round(coef(m)[1], digits = 1)), 
                        b = as.numeric(round(abs(coef(m)[2]), digits = 1)), 
                        r2 = as.numeric(round(summary(m)$r.squared, digits = 2))))
  as.character(as.expression(eq))
}

getr2 <- function(y, x){
  m <- lm(y ~ x)
  r2 = format(summary(m)$r.squared, digits = 3)
}

df <- group_by(df, model) %>%
  mutate(label=lm_eqn(y, x))

p1 <- ggplot(df, aes(x=gdpcap_2010, y=(y/(2100-2010)), label=label)) + geom_point(size=0.5) +
  geom_smooth(method = "lm", se=FALSE, color="red", formula = y ~ x) +
  facet_wrap(~model, nrow=3, ncol=5) +
  geom_text(x=3.5, y=0.07, size=4, parse=TRUE) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_x_log10(labels = comma_format(scale=0.01)) + 
  theme_bw() + xlab("GDP/capita in 2010 (in hundreds)") + ylab("Average yearly growth rates (SSP, 2010-2100)") +
  theme(strip.text.x = element_text(size=14),
        axis.title = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.text.x = element_text(size=14, hjust=1, angle=25))

ggsave("results/figure2.tiff", plot=p1, width=7.5, height=6, dpi=300, compression="zip")